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Abstract. The Jeffcott equations are a system of coupled differential 


equations that represent the behavior of a rotating shaft. This is a simple 
model that allows investigation of the basic dynamic behavior of rotating 
machinery. Nonlinearities can be introduced by taking into consideration 
deadband, side force, and rubbing, among others. 

In this paper we study the properties of the solutions of the Jeffcott 
equations with deadband. In particular, we show how bounds for the solutions 
of these equations can be obtained from bounds for the solutions of the 
linearized equations. By studying the behavior of the Fourier transforms of 
the solutions, we are also able to. predict the onset of destructive 
vibrations. These conclusions are verified by means of numerical solutions 
of the equations, and of power spectrum density (PSD) plots. 

This study offers insight into a possible detection method to determine 
pump stability margins during flight and hot fire tests, and was motivated by 
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the need to explain a phenomenon observed in the development phase of the 
cryogenic pumps of the Space Shuttle, during hot fire ground testing; namely, 
the appearance of vibrations at frequencies that could not be accounted for 
by means of linear models* 
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1. Introduction 


H. H. Jeffcott [10] was one of the first to study the vibration 
characteristics of an unbalanced, uniform, flexible shaft supported by 
bearings. He did so by considering a linear system of differential equations 
of the form 


y" + Cy' + Ay ■ F cos wt 


( 1 ) 


z" + Cz' + Az = F sin wt 

where the differentiation Is with respect to the parameter t, and the shaft 
is assumed to rotate along the x-axls with angular velocity w, y and z 
describe the displacement of the center of the shaft, and the coefficients 
have the following physical interpretation: C ■ C^/m, A ” K^/m, K * K^/m, 

2 

and F » uw , whre m is the mass of the shaft, C is the seal damping, 

s 

K and K. are the seal and bearing stiffnesses, and u is the displace- 
8 D 

ment of the shaft's center of mass from the geometric center. 
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In [12] Yamamoto studied the effect of deadband (l.e. the clearing 


between housing and bearing), but his treatment was not rigorous. Other 
works dealing with nonl inearl ties Include [3], [4], [8] and [llj. In [7], Day 
used the method of multiple scales to gain new Insight into the properties of 
the solutions. He discovered a frequency, which he termed "nonlinear natural 
frequency” that appears in the PSD plots of solutions of the nonlinear model 
and Is absent from -the PSD plots of solutions of the linear model. The 
nonlinear natural frequency seems to have been observed during early ground 
testing of both LOX and fuel pumps of the second stage Main Engine of the 


Spac.e Shuttle but, .until now, there has been no explanation of its origin.. 

If r - (y 2 + z 2 ) 1/2 , 6 is the deadband, K - K b /®» where K b Is the bearing 
stiffness, B * Q /m, with Q denoting the cross-coupling stiffness of the 

8 S 



then the model studied by Day can be described by the system 


y" + Cy' + [A + K(l-h) ]y + Bz - fj(t), 
z" + Cz' - By + [A + K(l-h) ] z - f 2 (t). 


XXXVI-2 



where f^(t) ** F cos wt, f^(t) m F sin wt, and K(l—h) is the nonlinearity 
associated with the deadband (see Fig, 1), Note that (1) is a particular 
case of (3). 

In this paper we shall explain the nature of the nonlinear natural 
frequency and apply our conclusions to the signature analysis of the 
nonlinear Jeffcott model described by (3), where f^(t) and f^(t) are 
arbitrary bounded and continuous functions, B, C, K, and 6 are positive, 
and A and t are nonnegative. 

2. Properties of the solutions of Jeffcott' s equations. 

2.1 Existence, uniqueness, and a representation formula. 

If Xj - y, x 2 - y’, x 3 - z, x 4 - z', r - (x 2 + x 2 ) 1/2 and h(t) Is 

given by (2), then (3) is equivalent to 
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I 


X- - - [A + K(l-b)lx, - Cx_ - Bx_ + f.(t) 

2 1 Z J 1 (4) 

X 3 “ X 4 

x^ - - Bx x - [A + K(l-n)]x 3 - Cx^ + f 2 <t) 
or, more concisely, 

I 

x - £<X, t). 

Since ^(x,t) is continuous, and satisfies a Llpschltz condition on x, from 
standard existence and uniqueness theorems ( cf . , eg* [5]), we know that 
every initial value problem for (4) has a unique solution. Thus,, we also 
infer that every initial value problem for (3) has a unique solution. Let 
v - y + iz, f(t) - f 1 (t)+if 2 (t), and M - A+K - IB; then (3) is also 
equivalent to 

v” + Cv' + Mv - Khv - f(t). 

Before studying (5), let us first consider a linear system of the form 


v” + Cv* + Mv = g(t). 


( 6 ) 
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Then 


v “ C i + C 2 exp(X 2 t) + v p , where C 1 exp(X, t)+C 2 exp( A t) is a 


solution of 


v" + Gv' + Mv - 0, 


( 7 ) 


and therefore A^ - (1/2) l- C + (C* - 4M) 1/2 ]. If Q - c 2 - 4(A + K), 


straightforward computation shows that A 


1 * a + 18, A » a' - ig, where 


a"l/2 


8 ■ 8 [-Q + (Q* + 16 B 2 ) 1 ^ 2 ] 1 ^ 2 , 


,-l 


a - [6 * B - C]/2 , a' - -[g _1 B + C]/2, 


( 8 ) 

(9) 


and therefore 


a -' - -g” 1 B + a 


( 10 ) 


Applying, e.g. [5, Theorem 6.4] we readily deduce that the Green's function 


of the differential operator (d 2 /dt 2 ) + C(d/dt) + M is G(t - s), where 


G(t) - (A 2 - Aj) 1 [exp( A^ t) + exp(A 2 t)], i.e. / G(t-s)q(s)ds is a 


particular solution of (6), and therefore 
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If in particular g(s) = K h(s) v(s) + f(s), then (6) reduces to (5), and we 


have: 

v “ c, exp(A.t) + c 0 exp (At) + / G(t~s)f(s)ds + Kj G((t—8)h(s)v(8)ds« 
1 1 2 2 0 0 

In other words, 

v(t) = u(t) + P(t), 

where 

u(t) ■ c. exp(A.t) + c„ exp(A t) + f t G(t-s)f(s)ds 

i 1 z z 0 

is a solution of the linear differential equation v" + Cv' + Mv - f(t), 
(henceforth called the linear part of (5)), and 
p(t) « K/q G(t-8) h(s)v(s)ds. 

Note that we have a closed form formula for u, whereas P(t) is 
expressed in terms of the unknown function v(t). Our analyses will be 
on a study of the properties of the perturbation term P(t). 


(ID 


( 12 ) 


(13) 


based 
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2.2 Bounds 


The definition (2) of h(s) implies that |h(t)v(t)| < 6; thus, 


|p(t,| < K6 / |G( t-s) | ds _< 


K6 (B 2 B 2 + 4 8 2 ) *^ 2 / [exp (a(t-s)) + exp(a' (t-s))]ds. 

0 


- 1 , 


Since (9) and (10) imply that if a - 0, then B ■ B/C, and when a » 8 B 


then a' ■ 0, we have 


|P(i 


t) < 


^K6(B~ 2 B 2 + 4B 2 ) 1 ^ 2 [(l/a)(exp(at)-l) + (1/a' )(exp(a' t)-l) ] 

if a * 0, 6~ 1 B 

9 

K6(C 2 + 4B 2 /C 2 ) 1/2 [ t + (2/C) (l-exp[ (-C/2) t] ) (14) 

if a * 0, 

K6(b" 2 B 2 + 4B 2 )“ 1/2 [(B/B)(exp[(B/B)t] -1) + t] 


if a - B _1 B 


From (11), (12), (13), and (14) we derive the following conclusions: 
1. If a < 0 and G(t-s) f(s)ds| ^ M, then the steady state 


solution v of (5) satisfies the following inequality: 

OO 


|v»j + K6(B“ 2 B 2 + 4B 2 )“ 1/2 |l/a + 1/a* | 


XXXVI- 7 


2. If a « 0, the perturbation term P(t) can grow at most linearly. 

3. If a > 0, the order of growth of P(t) cannot exceed exp(at); note 
that the order of magnitude of all nonzero solutions of (7) cannot exceed 
exp(at). 

Thus, since we have assumed that f(t) is bounded, the study of the 
boundedness of the solutions of (5) reduces to the study of the boundedness 
of the solutions of its homogeneous part. If a < 0 we shall say that (5) 
is stable, if a > 0 that (5) is unstable, and if a - 0 that (5) has 
reached the stability boundary. This nomenclature is consistent with that 
used for linear systems (cf. [9, pp. 83, 84]). 

2.3 Estimates for 8. 

1/2 

In this section we will prove that 8 is between B/C and (A+K) 

The importance of this observation will become clear in the sequel. In what 
follows, let y - A + K and t, - (B/C)^; thus Q * - 4y. 
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Assume that a < 0; then from (9) we see that B/C < 6. Squaring and 
applying (8), we have 

5 < (1/8) [- Q + (Q 2 + 16 B 2 ) 1/2 ] . 

Thus, £ satisfies the inequality 4£ 2 + Q5 - B 2 < 0, which can be written as 

45 2 + [C 2 - 4y] £ - B 2 < 0, or 

C - 4y < (B - 4? ) K m B 2 £ ^ - 4 £. Since £ =* (B/G) 2 , we have that 

2 2 

C - 4y < C - 45, and therefore Y > 5, i.e. 

2 2 
o < C y 

Thus, 

16y 2 + C 4 - 8C 2 Y + 16B 2 < 16y 2 + C 4 + 8C 2 y, 

1 » 6 « 

Q 2 + 16B 2 < (4y + C 2 ) 2 , 

and therefore 

[Q 2 + 16 B 2 ] 1/2 < 4y + C 2 . 
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Subtracting Q from both sides of this inequality we see that 
86 2 ■ - Q + [Q 2 + 16B 2 ] 1/2 < 8y, 

1/2 

and we conclude that 6 < (A + K) • 

The cases a - 0 and a > 0 are treated similarly, and we shall omit the 

details. The conclusions are the following: 

1/2 

1. If a < 0, then B/C < 8 < (A + K) , and a' < 0. 

1/2 

2. If a - 0, then B/C » 8 - (A + K) , and a' < 0. 

1 /2 

3. If a > 0, then (A + K) A/ < 8 < B/C. 

From these conclusions we also Infer that If f(t) is bounded, then (5) is 

1/2 

stable- if and only If B/C < (A + K) . 

2.4 Resonance. 

From the results of 2.3 it is clear that (5) is in resonance if and only if 
its linear part is in resonance. If for example f(t) ■ Fq exp(iwt), then we 
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readily see that the linear part of (5) has a particular solution of the form 


Aq exp(iwt), where 

A Q - F q / [ (A + K - W 2 ) + i(w-e)]. 

Since the denominator in the preceding formula vanishes if and only if 

1/2 1/2 

8 * w and w ■ (A + K) , we deduce that 8 * (A + K) , and therefore 
that a - 0. Thus (5) can be in resonance only on the stability boundary. 

3. Harmonic Analysis of the solutions. 

3.1 Introduction 

In practice, the coefficients of (5) and, in general, the equations that 
describe the movement of rotating machinery, are imperfectly known. The 
approach taken is to sample the system response over a time interval (in our 
case, that would mean measuring y(t^) and zCt^), i ■ 0, ..., N, where the 
t^ are equally spaced poins), and to approximate the Fourier transforms of 
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y(t) and z(t) by means of the Discrete Fourier Transforms of the sequences 


{yCt^)} and {z(t 1 )}. (See, e.g. [1]). The absolute values of the 
coefficients In the Discrete Fourier expansions are then plotted on graphs 
called Power Spectrum Density (PSD) plots, which represent the response of 
the mechanical system at different frequencies. One then tries to determine 
the condition of the mechanical system by an examination of these plots. 

This Is known as "signature analysis". (See, e.g. Collacott [6].) In this 
section we examine the properties of the Fourier transforms of the solutions 
of (5), whereas in section A we show, by means of examples, how to apply 
these conclusions to the signature analysis of the system. 

From now on, we shall assume that a < 0. 

3.2 Properties of the continuous Fourier transform. 

Let G(t), v(t), u(t), and P(t) be defined to equal 0 for t < 0. Then 
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(11) Is valid on (-<*>, ”). If 

q(t) - h(t) v(t), 

it is dear that q(t) vanishes for t < 0. Thus, from (13) we readily see 
that 

00 

P(t) - K / G(t-x) q(x)dx - K(G * q)(t), (13) 

— oo 

where "**' denotes the convolution product. If p^(t) » exp(X^t), 

P 2 (t) * exp(X 2 t) for t ^ 0, and equal zero for t < 0, and u p (t) denotes 
a particular solution of the linear part of (5), then (11) can be written in 
the form 

v(t) - c x p A (t) + c 2 P 2 (t) + u p (t) + P(t). (16) 

Note, moreover, tha t 

G( t) » (X 2 - X 1 ) -1 [p 1 (t) + p 2 (t)] and therefore 

P(t) - (X 2 - Xj)" 1 [(p x + p 2 ) * q] (t). (17) 

Let F denote the Fourier transform operator; thus, if g(t) is 
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00 


lntegrable on (-°°, °°)» then F[g](s) ■ / g(t) exp(-ist)dt. 

— oo 

In particular, when a < 0 we have 

Fl Pl )(s) - l/[i(s - B) - a] U8) 

and 

F[p 2 ](s) - 1/ [ i( 8 + B) - a + 6 _1 B] . (19) 

If u (t) and q( t) are lntegrable on (- 00 , “), from (16) and (17) we 
P 

see that 

F[v] - Cj_ FlpJ + c 2 F[p 2 ] + Flu p ] + (X 2 - X 1 )" 1 (F[p 1 ] + Flp 2 J) F(q). 

Since Ftpj^KB) diverges as a ♦ o”, we therefore conclude that if a 
approaches 0 from the left, then the graphs of the real and imaginary parts 
of F[v] will exhibit increasingly large spikes at 6, where a and B are 
linked by (9). At first glance, this does not appear to be very useful, 
since In most applications u p will not be lntegrable on (— , 00 ) (as for 

example when f(t) - F exp(lwt)). We shall now show that the range of 
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validity of our conclusions can be greatly extended If we consider windows 


3.3 Windowing and the nonlinear natural frequency 
3.3.1 Analysis of the transient terras 

In practice, Fourier transforms are computed for samples taken over a time 
interval of the form (a,b), (called a "window”), where a is in general 
larger than zero. Let g^ a,b \t) ■ g(t) if a _< t _< b, let g^ a * b \t) 
equal zero otherwise, and let g^ - g^* b ^. Clearly 
F [p(a,b)j( a ) _ [exp(Aj - si)b - exp(A 1 - sDa]/^ - si), 
and 

F[p 2 a,b ^ ] (s) - [exp(X 2 - 8i)b - exp(X 2 - si)a]/(A 2 - si). 

Since Aj - 8i - o, we see that F[pj a,b ^](8) diverges as a + 0 ; thus, 

also the graphs of the real and imaginary parts of F[v^ a,b ^] should have 

spikes at 8. However, since lim F[pj a,b ^](s) ■ lim F[p^ a,b ^](s) « 0, we 

a+» a+® 
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conclude that for fixed a and sufficiently large a, these spikes may be 


detected only If a Is extremely close to 0. (Whether they will be 
detected at all depends on the numerical stability of the computations). 

Thus, in order to obtain useful data we have to analyze the Fourier transform 
of the perturbation term P(t). 

3.3.2 Analysis of the perturbation term 

qo £ 

Let P (t) ■ K / G(t-x) q^(x) dx « / G(t - x) q^ b \x) dx. From (15) we 
b 0 

see that if t <_ b, then P b (t) “ P(t), whereas for t ^ b, 

P. (t) ■ f b G(t - x) q(x) dx. Thus, 
b 0 

F[P (b) ](s) - F[P b ](s) - K I b (s), (20) 

where 

oo 5 

I (s) * / exp(-sti) / G(t - x) q(x) dx dt. 

D b 0 

We can write 1. in the form 
b 

h “ (x 2 " '1 )_1 lI(A i } ■ 

with 
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KX) - /” exp(-sti) /J exp [ A ( t - x)] q(x) dx dt. 

If Re(A) < 0, reversing the order of integration we have: 

1(A) ■ /q exp(-Ax)q(x)dx J exp[(A - si)]dt 

= M(b, A) exp [ ( A - si)b]/(A - si), 
b 

where M(b, A) = / exp(-Ax)q(x) dx. 

0 

Thus, 

“ (A 2 - Aj) 1 M(b, A 1 ) exp [ ( A ^ - si)b]/(A 1 - si) 

( 21 ) 

+ (A 2 - Aj) 1 M(b, A 2 ) exp[A 2 - si)b]/(A 2 - si). 

Moreover, since q^(x) is .of bounded support it is integrable. Thus, 
since p b (t) * K (G * q^)(t), we know that 

F[P b Ks) « K F[G](s) F[q (b) J(s) (22) 

Since G(t) ■ (A 2 - A^) 1 tp^(t) + P 2 (t)], from (18) and (19) we have: 
FtG](s) - (A 2 - A^" 1 [(A x - si)" 1 + (A 2 - si)" 1 ]. 
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Combining (20), (21), and (22), we thus obtain: 

F[P (b) ](s) - 

K(X 2 - Xj) -1 (F[q (b) l(s) - M(b, X^expt^ - si)b)}/(X 1 - si) + 

K(X 2 - X^" 1 (F[q^ b ^](s) - M(b, X 2 ) expl(X 2 - si)b]}/(X 2 - si). 

Thus, since F[P (a » b) ](s) - FlP (b) ](s) - F[P (a) ](s), 
setting 

Q(a, », X, s) - (23) 

F[q( a » b )]( 8 ) _ M(b, X)exp[ (X - si)b] + M(a, X)exp[(X - si)a], 

we conclude that 

p[p( a * b ) ] ( 8 * ■ Q(a, b, Xj, s)/(X 1 - si) + Q(a, b, X 2 , s)/(X 2 - si). (24) 
If Q(a, b, X x , e) * 0, we conclude that F[P (a,b> l(8) will diverge as 
a 0 . 

We shall now show that, for any e > 0, the functions Q(a, b, X^, s) 
and Q(a, b, X 2 , s) are bounded, uniformly on s, provided that 
- e < a < 0. From (23) it is clear that 
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Now, the 


|q(a, b, X, s)| < | F[q (a,b) ] (s) | + |M(b, X)| + |M(a, X)|. 
definition of G(t) implies that |G(t)| — ” ^l| Since 

|h(t)v(t) |<6 , It Is clear from (13) that |P(t)| 2K|X 2 ” ^1 1 

Assuming that |f(t)| S. M i * * 8 readily seen from (11) and (12), that 

| v( t) | < JCj | + |c 2 | + 2 |x 2 - X 1 |“ 1 (kfi + Mj) t. 

The constants and C 2 depend on the Initial conditions. We shall now 

show that, for the same set of Initial conditions and any t > 0, and C 2 
are bounded, uniformly on a. 

Assume that v(0) ■ Vq, and v'(0) ■ v^. From (11), (12), and (13), it 
is clear that 

Cj + C 2 - , 0 . (25) 

On the other hand, since 

t t 

(d/dt) / G(t-x)g(x)dx « G(0)g( t) + / (3/9t)G(t-x)g(x)dx, 

0 0 

(cf., e.g. Bar tie [2]), differentiating (11) we have 


XXXVI- 19 


X 1 C 1 + X 2 C 2 + g ^°^ 8(0) * Vj, where g(t) ■ f(t) + Kh(t)v(t), and from (25) 


we obtain 


X 1 C 1 + X 2^ V 0 ” C P + ■ Vj. 

If |f(t)| S. M i » an< * M - Mj + K6 vq, it is easy to see that |g(0)| M. 
Thus, since G(0) ■ (X^ - X^) X , we deduce that 



Since we are assuming that B and C stay positive, from (9) and- (10) we 
see that, as a 0 , X^, X^ and (X^ - ^) remain bounded. Thus, also 
c^ and c 2 remain bounded, and from (24) we see that if K reamlns bounded, 
then for every e > 0 there are constants A Q , B^, that do not depend on a, 
such that |v(t)| <. A Q + B Q t for a in (-e, 0). Applying this Inequality 
it is now easy to see that for any c > 0, Q(a, b, X^, s) and 
Q(a, b, X s) are bounded, uniformly on a, provided that -c < a < 0. 

Thus, we conclude that there are constants * K^(e), and K 2 (e), such that 

|F[P (a,b) ](s)| < K x / |X X - s i | + K 2 /|X 2 - s?|. 

This means that the only value for which F[P^ a,b ^](s) diverges as a * 0 
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is s ■ 8. Also, there is no obvious reason why F[P^ a,b ^](8) should vanish 
as a » (provided that we keep the difference b - a constant). 

3.3.3 Conclusions 

In summary, we have shown that F[v^ a,b ^](8) diverges as o + O - , that 
8 is the only value for which this may happen, that for a negative and 
constant, but sufficiently close to zero, the graphs of the absolute values 
of the real and imaginary parts of F[v^ a,b ^](s) will have spikes at s ■ 8, 

and that the magnitude of these spikes need not decrease with time (i.e., as 

\ * 
a + 00 ) • 

In [7, p. 784], Day equates the nonlinear natural frequency with the 
ratio Q /C (which, In our notation, equals B/C). Later on, (on p. 786), 

S 8 

he notes that the nonlinear natural frequency is actually not B/C, but a 
number close to it. In this pape** we have gone one step further and shown 
that 8 (i.e. the transient frequency of the linear part of (5)), and the 

nonlinear natural frequency are one and the same. This is a very surprising 
result. 
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On the practical side, our conclusions suggest that, all other things 


being equal, the introduction of nonlinearlties (as Induced, e.g. by 
deadbands) In a mechanical system, may give an earlier warning of the 
approach to the instability boundary. 

4. Examples 

We now study the behavior of the solutions of (5) for C ■ 240, A * 0, 

6 - 0.0000285, K - 1,305,000, f(t) - uw 2 exp(iwt), u ■ 0.00006915, and 
w ” l,000trs, where s will vary. We also make the realistic assumption 

f 

-based on empirical data- that the bearing stiffness changes with the forcing 
frequency w, by setting B - 60w. Let f (the "critical frequency") be 

defined to equal (A-HC) 1 ^ 2 /(2 ti) - K 1 ^ 2 /(2ir). We readily see that 

f A 181.8 Hz., and that the value of s that corresponds to 
c * 

f is s A 1.4545. For s » 10/3 our example reduces to Example 1 of [7]. 
c c 
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Let fj ■ w/2tt denote the forcing frequency; clearly f^ ■ 500 s. We know 
that (5) will be stable for s < s c , and unstable for s > s c> In Figs. 2 
through 5 we show, for various values of s, plots of the numerical solutions 
of (5), (obtained by a fourth order Runge-Kutta algorithm), and of PSD's for 
the real part y and imaginary part z of v. The solution plots are for 
0.1 < t < 0.256, and the PSD plots for the window [0, 0.256]. Note that all 
the PSD plots have two distinct spikes: one corresponding to the forcing 
frequency, and one corresponding to the nonlinear natural frequency. For s 
- 0.5 (i.e. far from the stability boundary), the forcing function f(t) 
dominates tlv perturbation term P(t).. Thus, the solution is nearly circular 
(we see a thick circular curve; the thickness is caused by P(t)), and the 
PSD plots exhibit larger spikes for the forcing frequency than for the 
nonlinear natural frequency. As s Increases, the solution becomes annular, 
and the nonlinear natural frequency begins to dominate. Finally, for s > s^ 
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the solutions begin to diverge. Since the PSD plots are obtained by 
approximating Fourier transforms by discrete Fourier transforms, (which are 
intrinsically bounded) they show no obvious qualitative difference when 
compared with plots for values of a close to, but smaller than, s c< 

Figures 6 and 7 snow PSD plots for s * 1.2 and s - 1.3 and various 
windows. Note that if we compare the PSD plots for 0 < t < 0.256 (Figs. 3 
and 4) and 0.256 < t < 0.512 (Figs. 6 and 7), we see a large decrease in the 
height of the spike that corresponds to the nonlinear natural frequency, and 
a very small decrease when we compare the plots for 0.256 < t < 0.512 and 
l;D24 < t < 1.28, but there is no change, in magnitude on subsequent windows. 
This is due to the disappearance of the transient terms F [ p ^ ] and F^j* 
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